Clustering Statistics
C. Savonen CCDL for ALSF
12/17/2018
Purpose: compare normalization methods for Darmanis et al data
# Load in functions
source(file.path("scripts", "CZI_functions.R"))
# These are the necessary packages
packages <- c("Rtsne", "NMI", "caret")
# Check if these packages are installed and install them if they aren't
lapply(packages, function(package) {
if (!(package %in% installed.packages())) {
install.packages(package)
}
})## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULLImport data that was previously obtained from submitting counts that were formatted by arnon_asap_data_prep.Rmd to ASAP
# Get the file names of all the normalized files
normalized.files <- dir("darmanis_data/normalized_darmanis")
# Read in each of the normalization files
normalized.data <- lapply(normalized.files, function(file) {
readr::read_tsv(file.path("darmanis_data/normalized_darmanis",
file))
})## Parsed with column specification:
## cols(
## .default = col_double(),
## gene = col_character()
## )## See spec(...) for full column specifications.## Parsed with column specification:
## cols(
## .default = col_double(),
## Genes = col_character()
## )## See spec(...) for full column specifications.## Parsed with column specification:
## cols(
## .default = col_double(),
## Genes = col_character()
## )## See spec(...) for full column specifications.## Parsed with column specification:
## cols(
## .default = col_double(),
## Genes = col_character()
## )## See spec(...) for full column specifications.## Parsed with column specification:
## cols(
## .default = col_double(),
## Genes = col_character()
## )## See spec(...) for full column specifications.## Parsed with column specification:
## cols(
## .default = col_double(),
## Genes = col_character()
## )## See spec(...) for full column specifications.# Keep all the gene lists
genes <- lapply(normalized.data, function(x) x[,1])
# Keep the names with it.
names(normalized.data) <- gsub("\\.tab|darmanis_", "", normalized.files)# Read in the data
meta <- readr::read_tsv(file.path("darmanis_data", "GSE84465_meta.tsv"))## Parsed with column specification:
## cols(
## .default = col_character(),
## channel_count = col_double(),
## taxid_ch1 = col_double(),
## contact_phone = col_double(),
## `contact_zip/postal_code` = col_double(),
## data_row_count = col_double(),
## `plate id:ch1` = col_double(),
## `tsne cluster:ch1` = col_double()
## )## See spec(...) for full column specifications.Retrieve cell type and plate batch for these samples:
# Get sample names from columns
samples <- colnames(normalized.data[[1]])[-1]
# Keep metadata only for the samples we have
meta <- meta[match(samples, meta$geo_accession), ]
# Extra cell types info as it's own vector
cell.types <- as.factor(meta$`cell type:ch1`)
plate.batch <- as.factor(gsub("plate id: ", "", meta$`characteristics_ch1.1`))Run tsne clustering for each dataset:
# Run tsne on each dataset and extract
tsne <- lapply(normalized.data, function(dat) {
tsne.res <- Rtsne::Rtsne(t(dat[,-1]),
check_duplicates = FALSE)
tsne.res <- tsne.res$Y
names(tsne.res) <- samples
return(tsne.res)
})Plot all the tsne’s for all the datasets with labeling with batch info as well as cell type info.
# Make a function for plotting tsne's and using
tsne.plot <- function(dat, var, name = "name"){
# This function is to plot tsne data and label it by a variable
# Args:
# dat: a data.frame with two columns of data
# var: vector contains metadata labels
# name: a character string for the plot title and for to save as png
# Returns:
# png and plot in Rmd of the provided data with labels of the metadata provided
# convert celltype into numbers
colz <- colors(distinct = TRUE)[runif(length(levels(var)), min = 1,
max = length(colors(distinct = TRUE)))]
plot(dat, pch = 21, bg = colz[var], main = name, ylab = "tsne dim 2",
xlab = "tsne dim 1");
legend(x = "bottomleft", legend = levels(var), fill = colz, cex = 0.6)
}Plot all the tsnes
# Plot em with cell type and plate batch labels
lapply(tsne, function(dataset) {
# Get normalization method
set.name <- names(tsne)[parent.frame()$i[]]
# Make plots for cell type and plate batch
cell.type.plot <- tsne.plot(dataset, cell.types, name = set.name)
plate.batch.plot <- tsne.plot(dataset, plate.batch, name = set.name)
# Save the plots to pngs
png(paste0("tsne_", set.name, "cell_type.png"))
cell.type.plot
dev.off()
png(paste0("tsne_", set.name, "plate_batch.png"))
plate.batch.plot
dev.off()
# Print out the plots in the Rmd
cell.type.plot
plate.batch.plot
})
## $counts
## $counts$rect
## $counts$rect$w
## [1] 7.378351
##
## $counts$rect$h
## [1] 101.8396
##
## $counts$rect$left
## [1] -23.18612
##
## $counts$rect$top
## [1] 77.82416
##
##
## $counts$text
## $counts$text$x
## [1] -21.00673 -21.00673 -21.00673 -21.00673 -21.00673 -21.00673 -21.00673
## [8] -21.00673 -21.00673 -21.00673 -21.00673 -21.00673 -21.00673 -21.00673
## [15] -21.00673 -21.00673 -21.00673 -21.00673 -21.00673 -21.00673 -21.00673
## [22] -21.00673 -21.00673 -21.00673 -21.00673 -21.00673 -21.00673 -21.00673
## [29] -21.00673 -21.00673 -21.00673 -21.00673 -21.00673 -21.00673 -21.00673
## [36] -21.00673 -21.00673 -21.00673 -21.00673 -21.00673 -21.00673 -21.00673
## [43] -21.00673 -21.00673 -21.00673 -21.00673 -21.00673 -21.00673 -21.00673
## [50] -21.00673
##
## $counts$text$y
## [1] 75.82730776 73.83045317 71.83359858 69.83674400 67.83988941
## [6] 65.84303482 63.84618024 61.84932565 59.85247106 57.85561648
## [11] 55.85876189 53.86190730 51.86505272 49.86819813 47.87134354
## [16] 45.87448896 43.87763437 41.88077978 39.88392520 37.88707061
## [21] 35.89021602 33.89336144 31.89650685 29.89965226 27.90279768
## [26] 25.90594309 23.90908850 21.91223392 19.91537933 17.91852474
## [31] 15.92167015 13.92481557 11.92796098 9.93110639 7.93425181
## [36] 5.93739722 3.94054263 1.94368805 -0.05316654 -2.05002113
## [41] -4.04687571 -6.04373030 -8.04058489 -10.03743947 -12.03429406
## [46] -14.03114865 -16.02800323 -18.02485782 -20.02171241 -22.01856699
##
##
##
## $DESeq2
## $DESeq2$rect
## $DESeq2$rect$w
## [1] 10.77448
##
## $DESeq2$rect$h
## [1] 77.43283
##
## $DESeq2$rect$left
## [1] -40.08979
##
## $DESeq2$rect$top
## [1] 59.45075
##
##
## $DESeq2$text
## $DESeq2$text$x
## [1] -36.90728 -36.90728 -36.90728 -36.90728 -36.90728 -36.90728 -36.90728
## [8] -36.90728 -36.90728 -36.90728 -36.90728 -36.90728 -36.90728 -36.90728
## [15] -36.90728 -36.90728 -36.90728 -36.90728 -36.90728 -36.90728 -36.90728
## [22] -36.90728 -36.90728 -36.90728 -36.90728 -36.90728 -36.90728 -36.90728
## [29] -36.90728 -36.90728 -36.90728 -36.90728 -36.90728 -36.90728 -36.90728
## [36] -36.90728 -36.90728 -36.90728 -36.90728 -36.90728 -36.90728 -36.90728
## [43] -36.90728 -36.90728 -36.90728 -36.90728 -36.90728 -36.90728 -36.90728
## [50] -36.90728
##
## $DESeq2$text$y
## [1] 57.9324555 56.4141647 54.8958738 53.3775830 51.8592921
## [6] 50.3410013 48.8227105 47.3044196 45.7861288 44.2678379
## [11] 42.7495471 41.2312563 39.7129654 38.1946746 36.6763837
## [16] 35.1580929 33.6398021 32.1215112 30.6032204 29.0849295
## [21] 27.5666387 26.0483479 24.5300570 23.0117662 21.4934753
## [26] 19.9751845 18.4568936 16.9386028 15.4203120 13.9020211
## [31] 12.3837303 10.8654394 9.3471486 7.8288578 6.3105669
## [36] 4.7922761 3.2739852 1.7556944 0.2374036 -1.2808873
## [41] -2.7991781 -4.3174690 -5.8357598 -7.3540506 -8.8723415
## [46] -10.3906323 -11.9089232 -13.4272140 -14.9455049 -16.4637957
##
##
##
## $log2
## $log2$rect
## $log2$rect$w
## [1] 7.881444
##
## $log2$rect$h
## [1] 86.66107
##
## $log2$rect$left
## [1] -27.75925
##
## $log2$rect$top
## [1] 64.34185
##
##
## $log2$text
## $log2$text$x
## [1] -25.43127 -25.43127 -25.43127 -25.43127 -25.43127 -25.43127 -25.43127
## [8] -25.43127 -25.43127 -25.43127 -25.43127 -25.43127 -25.43127 -25.43127
## [15] -25.43127 -25.43127 -25.43127 -25.43127 -25.43127 -25.43127 -25.43127
## [22] -25.43127 -25.43127 -25.43127 -25.43127 -25.43127 -25.43127 -25.43127
## [29] -25.43127 -25.43127 -25.43127 -25.43127 -25.43127 -25.43127 -25.43127
## [36] -25.43127 -25.43127 -25.43127 -25.43127 -25.43127 -25.43127 -25.43127
## [43] -25.43127 -25.43127 -25.43127 -25.43127 -25.43127 -25.43127 -25.43127
## [50] -25.43127
##
## $log2$text$y
## [1] 62.6426109 60.9433742 59.2441375 57.5449008 55.8456641
## [6] 54.1464274 52.4471906 50.7479539 49.0487172 47.3494805
## [11] 45.6502438 43.9510070 42.2517703 40.5525336 38.8532969
## [16] 37.1540602 35.4548235 33.7555867 32.0563500 30.3571133
## [21] 28.6578766 26.9586399 25.2594032 23.5601664 21.8609297
## [26] 20.1616930 18.4624563 16.7632196 15.0639828 13.3647461
## [31] 11.6655094 9.9662727 8.2670360 6.5677993 4.8685625
## [36] 3.1693258 1.4700891 -0.2291476 -1.9283843 -3.6276211
## [41] -5.3268578 -7.0260945 -8.7253312 -10.4245679 -12.1238046
## [46] -13.8230414 -15.5222781 -17.2215148 -18.9207515 -20.6199882
##
##
##
## $scVLM
## $scVLM$rect
## $scVLM$rect$w
## [1] 10.90754
##
## $scVLM$rect$h
## [1] 75.06794
##
## $scVLM$rect$left
## [1] -33.20834
##
## $scVLM$rect$top
## [1] 54.67521
##
##
## $scVLM$text
## $scVLM$text$x
## [1] -29.98652 -29.98652 -29.98652 -29.98652 -29.98652 -29.98652 -29.98652
## [8] -29.98652 -29.98652 -29.98652 -29.98652 -29.98652 -29.98652 -29.98652
## [15] -29.98652 -29.98652 -29.98652 -29.98652 -29.98652 -29.98652 -29.98652
## [22] -29.98652 -29.98652 -29.98652 -29.98652 -29.98652 -29.98652 -29.98652
## [29] -29.98652 -29.98652 -29.98652 -29.98652 -29.98652 -29.98652 -29.98652
## [36] -29.98652 -29.98652 -29.98652 -29.98652 -29.98652 -29.98652 -29.98652
## [43] -29.98652 -29.98652 -29.98652 -29.98652 -29.98652 -29.98652 -29.98652
## [50] -29.98652
##
## $scVLM$text$y
## [1] 53.2032888 51.7313685 50.2594482 48.7875279 47.3156075
## [6] 45.8436872 44.3717669 42.8998466 41.4279263 39.9560060
## [11] 38.4840857 37.0121654 35.5402451 34.0683248 32.5964045
## [16] 31.1244842 29.6525639 28.1806436 26.7087233 25.2368029
## [21] 23.7648826 22.2929623 20.8210420 19.3491217 17.8772014
## [26] 16.4052811 14.9333608 13.4614405 11.9895202 10.5175999
## [31] 9.0456796 7.5737593 6.1018390 4.6299186 3.1579983
## [36] 1.6860780 0.2141577 -1.2577626 -2.7296829 -4.2016032
## [41] -5.6735235 -7.1454438 -8.6173641 -10.0892844 -11.5612047
## [46] -13.0331250 -14.5050453 -15.9769657 -17.4488860 -18.9208063
##
##
##
## $TMM
## $TMM$rect
## $TMM$rect$w
## [1] 8.592886
##
## $TMM$rect$h
## [1] 80.45252
##
## $TMM$rect$left
## [1] -28.76301
##
## $TMM$rect$top
## [1] 57.58385
##
##
## $TMM$text
## $TMM$text$x
## [1] -26.22488 -26.22488 -26.22488 -26.22488 -26.22488 -26.22488 -26.22488
## [8] -26.22488 -26.22488 -26.22488 -26.22488 -26.22488 -26.22488 -26.22488
## [15] -26.22488 -26.22488 -26.22488 -26.22488 -26.22488 -26.22488 -26.22488
## [22] -26.22488 -26.22488 -26.22488 -26.22488 -26.22488 -26.22488 -26.22488
## [29] -26.22488 -26.22488 -26.22488 -26.22488 -26.22488 -26.22488 -26.22488
## [36] -26.22488 -26.22488 -26.22488 -26.22488 -26.22488 -26.22488 -26.22488
## [43] -26.22488 -26.22488 -26.22488 -26.22488 -26.22488 -26.22488 -26.22488
## [50] -26.22488
##
## $TMM$text$y
## [1] 56.0063494 54.4288490 52.8513485 51.2738480 49.6963476
## [6] 48.1188471 46.5413466 44.9638462 43.3863457 41.8088452
## [11] 40.2313448 38.6538443 37.0763439 35.4988434 33.9213429
## [16] 32.3438425 30.7663420 29.1888415 27.6113411 26.0338406
## [21] 24.4563401 22.8788397 21.3013392 19.7238387 18.1463383
## [26] 16.5688378 14.9913373 13.4138369 11.8363364 10.2588360
## [31] 8.6813355 7.1038350 5.5263346 3.9488341 2.3713336
## [36] 0.7938332 -0.7836673 -2.3611678 -3.9386682 -5.5161687
## [41] -7.0936692 -8.6711696 -10.2486701 -11.8261706 -13.4036710
## [46] -14.9811715 -16.5586719 -18.1361724 -19.7136729 -21.2911733
##
##
##
## $voom
## $voom$rect
## $voom$rect$w
## [1] 7.482998
##
## $voom$rect$h
## [1] 100.0519
##
## $voom$rect$left
## [1] -25.20722
##
## $voom$rect$top
## [1] 77.24933
##
##
## $voom$text
## $voom$text$x
## [1] -22.99693 -22.99693 -22.99693 -22.99693 -22.99693 -22.99693 -22.99693
## [8] -22.99693 -22.99693 -22.99693 -22.99693 -22.99693 -22.99693 -22.99693
## [15] -22.99693 -22.99693 -22.99693 -22.99693 -22.99693 -22.99693 -22.99693
## [22] -22.99693 -22.99693 -22.99693 -22.99693 -22.99693 -22.99693 -22.99693
## [29] -22.99693 -22.99693 -22.99693 -22.99693 -22.99693 -22.99693 -22.99693
## [36] -22.99693 -22.99693 -22.99693 -22.99693 -22.99693 -22.99693 -22.99693
## [43] -22.99693 -22.99693 -22.99693 -22.99693 -22.99693 -22.99693 -22.99693
## [50] -22.99693
##
## $voom$text$y
## [1] 75.2875328 73.3257312 71.3639296 69.4021280 67.4403264
## [6] 65.4785248 63.5167232 61.5549216 59.5931200 57.6313184
## [11] 55.6695168 53.7077151 51.7459135 49.7841119 47.8223103
## [16] 45.8605087 43.8987071 41.9369055 39.9751039 38.0133023
## [21] 36.0515007 34.0896991 32.1278975 30.1660958 28.2042942
## [26] 26.2424926 24.2806910 22.3188894 20.3570878 18.3952862
## [31] 16.4334846 14.4716830 12.5098814 10.5480798 8.5862782
## [36] 6.6244765 4.6626749 2.7008733 0.7390717 -1.2227299
## [41] -3.1845315 -5.1463331 -7.1081347 -9.0699363 -11.0317379
## [46] -12.9935395 -14.9553411 -16.9171428 -18.8789444 -20.8407460Make a function to do kmeans clustering and obtain Normalized Mutual Information “nmi”and Adjusted Rand Index “ari” statistics over iterations (default is 10).
kmeans_eval <- function(feature, metadata = metadata, iter = 10) {
# This function is used to perform iterative k-means clustering based on projected
# features for single cell data and then evaluate the performance according to
# Normalized mutual information (NMI) and adjusted rand index (ARI)
#
# Args:
# feature: a data.frame contains clustering dimensions
# metadata: vector contains cell type or other metadata information
# iter: number of interation for k-means clustering
# Returns:
# NMI and ARI results
#
# Convert celltype into numbers
metadata_num <- as.numeric(factor(metadata))
sample_id <- seq(1:nrow(feature))
# iterative k-means
nmi_score_all <- c()
ari_score_all <- c()
all_cluster <- list()
for(i in 1:iter){
# set k equal to the number of celltypes in the dataset
k <- length(unique(metadata))
# perform k means clustering
km <- kmeans(feature, k)
# true clusters
orignal_data <- data.frame(sample_id, metadata_num)
# predicted clusters
cl_data <- data.frame(sample_id, km$cluster)
# calculate NMI and ARI score
nmi_score <- NMI::NMI(orignal_data, cl_data)$value
ari_score <- mclust::adjustedRandIndex(km$cluster, metadata_num)
nmi_score_all <- c(nmi_score_all, nmi_score)
ari_score_all <- c(ari_score_all, ari_score)
}
# Compile all results into a data.frame
results <- data.frame(ari = ari_score_all, nmi = nmi_score_all)
return(results)
}Make a function to do K nearest neighbor clustering and obtain accuracy scores over iterations (default is 10)
knn_eval <- function(feature, metadata = metadata, iter = 10){
# This function performs knn based evaluation
# Args:
# feature: a data.frame contains clustering dimensions
# metadata: vector contains cell type or other metadata information
# iter: number of interation for knn clustering
# Returns:
# list of accuracy scores for each iteration of KNN
# Make the data into a data.frame:
feature <- data.frame("tsne" = feature, "metadata" = metadata)
# Split observations into groups
cv <- cvTools::cvFolds(nrow(feature), K = iter, R = 1)
# Create empty objects to store the performance information for each iteration
perf.eval <- list()
confusion.matrix <- 0
# Go through this iteration that k times
for (i in 1:iter) {
# Isolate samples for training the model
train <- feature[cv$subsets[-which(cv$which == i)], ]
# Isolate samples for testing the model
test <- feature[cv$subsets[which(cv$which == i)], ]
# Perform KNN model fitting
knn.fit <- caret::train(metadata~. , data = train, method = "knn",
trControl = caret::trainControl(method = "cv", number = 3),
preProcess = c("center", "scale"),
tuneLength = 10)
# Evaluate the model
knn.pred <- predict(knn.fit, newdata = subset(test, select = -c(metadata)))
perf.eval[[i]] <- round(cal_performance(knn.pred, test$metadata, 3), 2)
# Make the results into a matrix
matrix <- as.matrix(table(test$metadata, knn.pred, deparse.level = 0))
confusion.matrix <- matrix + confusion.matrix
}
# Get mean performance of cross validation
perf.eval <- dplyr::bind_rows(perf.eval)
accuracy <- perf.eval$accuracy
return(data.frame("knn" = accuracy))
}Run cell type analyses on all the tsne datasets
# Get knn and kmeans results for all tsne's of all datasets
cell.type.results <- lapply(tsne, function(dataset) {
knn.results <- knn_eval(dataset, metadata = cell.types)
kmeans.results <- kmeans_eval(dataset, metadata = cell.types)
data.frame(knn.results, kmeans.results)
})## Loading required package: lattice## Loading required package: ggplot2Run plate batch analyses on all the tsne datasets
# Get knn and kmeans results for all tsne's of all datasets
batch.results <- lapply(tsne, function(dataset) {
# knn.results <- knn_eval(dataset, metadata = plate.batch)
kmeans.results <- kmeans_eval(dataset, metadata = plate.batch)
# data.frame(knn.results, kmeans.results)
data.frame(kmeans.results)
})Set up a function for formatting the data into ggplot2 and subseuently plot the data
plot.results <- function(results.list, name = "results") {
# This function makes a boxplot of the cluster statistic results
# Args:
# results.list: a list of dataframes which each column contains a different
# statistic for 10 iterations
# name: name to use for the png to be saved and the plot title
# Returns: boxplots of the normalization method cluster statistics. Prints the
# plots and save the plot as a png
# Get meta info
meta <- names(unlist(results.list))
# Transform list into a dataframe:
ggplot.df <- data.frame("method" = stringr::word(meta, 1, sep = "\\."),
"test" = gsub("[0-9]*$", "",
stringr::word(meta, 2, sep = "\\.")),
"iter" = rep(1:nrow(results.list[[1]]),
length(results.list)*ncol(results.list[[1]])),
"values" = unlist(results.list))
# Make the plot
plot <- ggplot(data = ggplot.df, aes(x = method, y = values, fill = test)) +
geom_boxplot(position = position_dodge()) +
xlab("Normalization method") +
ggtitle(name) +
facet_wrap(~test)
# Save plot to png
ggsave(paste0(name, "_cluster_results.png"), width = 10)
# Print plot
plot
}Plot all the statistics on a boxplot
# Plot the cell type results
plot.results(cell.type.results, name = "Cell_type")## Saving 10 x 5 in image
# Plot the plate batch results
plot.results(batch.results, name = "Plate_batch")## Saving 10 x 5 in image
Session info:
sessionInfo()## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin17.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
## LAPACK: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libLAPACK.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] caret_6.0-81 ggplot2_3.1.0 lattice_0.20-35
##
## loaded via a namespace (and not attached):
## [1] tidyselect_0.2.5 purrr_0.2.5 reshape2_1.4.3
## [4] splines_3.5.1 colorspace_1.3-2 generics_0.0.2
## [7] stats4_3.5.1 htmltools_0.3.6 yaml_2.2.0
## [10] prodlim_2018.04.18 survival_2.42-6 rlang_0.3.0.1
## [13] e1071_1.7-0 ModelMetrics_1.2.2 pillar_1.3.0
## [16] NMI_2.0 withr_2.1.2 glue_1.3.0
## [19] bindrcpp_0.2.2 foreach_1.4.4 plyr_1.8.4
## [22] bindr_0.1.1 lava_1.6.4 robustbase_0.93-3
## [25] stringr_1.3.1 timeDate_3043.102 munsell_0.5.0
## [28] gtable_0.2.0 recipes_0.1.4 codetools_0.2-15
## [31] evaluate_0.12 labeling_0.3 knitr_1.20
## [34] class_7.3-14 DEoptimR_1.0-8 Rcpp_0.12.19
## [37] readr_1.3.0 scales_1.0.0 backports_1.1.2
## [40] ipred_0.9-8 hms_0.4.2 digest_0.6.18
## [43] stringi_1.2.4 Rtsne_0.15 dplyr_0.7.7
## [46] grid_3.5.1 rprojroot_1.3-2 tools_3.5.1
## [49] magrittr_1.5 lazyeval_0.2.1 tibble_1.4.2
## [52] crayon_1.3.4 cvTools_0.3.2 pkgconfig_2.0.2
## [55] MASS_7.3-51 Matrix_1.2-14 data.table_1.11.8
## [58] lubridate_1.7.4 gower_0.1.2 assertthat_0.2.0
## [61] rmarkdown_1.10 iterators_1.0.10 mclust_5.4.1
## [64] R6_2.3.0 rpart_4.1-13 nnet_7.3-12
## [67] nlme_3.1-137 compiler_3.5.1